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Abstract 

\ We derive a hierarchy of closures based on perturbations of well-known entropy-based 

i> closures; we therefore refer to them as perturbed entropy-based models. Our derivation 

£Lj reveals final equations containing an additional convective and diffusive term which are 

added to the flux term of the standard closure. We present numerical simulations for 
the simplest member of the hierarchy, the perturbed Mi or PM\ model, in one spatial di- 
mension. Simulations are performed using a Runge-Kutta discontinuous Galerkin method 
with special limiters that guarantee the realizability of the moment variables and the posi- 
tivity of the material temperature. Improvements to the standard Mi model are observed 
in cases where unphysical shocks develop in the Mi model. 

1 Introduction 

In this paper, we derive a new hierarchy of kinetic moment models in the context of frequency- 
integrated (grey) photon transport. These new models are perturbations of well known 
entropy-based models; we therefore refer to them as perturbed entropy-based or PEB models. 
We present numerical simulations for the simplest member of the hierarchy, the perturbed 
Mi or PM\ model, in one spatial dimension. In this setting, the PM\ model approximates 
the evolution of the photon radiation energy E and radiation flux F through a material 
medium with slab geometry. The photons interact with the material through scattering and 
emission/absorption processes. 

Entropy-based (EB) models have been studied extensively in areas such as extended ther- 
modynamics [201113], gas dynamics [Ml[Ml[321[33l[35l[38lE5], semiconductors [2H5H22EHI3H 
EE1I5I], quantum fluids [E1EI], radiation transport piIliniliaillEaBSESllSIllEllSlEaiSBl, 
and phonon transport in solids [21]. In the context of radiative transfer, entropy-based mod- 
els are commonly referred to as M/v models, where N is order of the expansion. The M\ 
model dates back to [43] . where it was first derived using Maxwell-Boltzmann statistics. For 
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problems with Bose-Einstein statistics, formal theoretical properties such as hyperbolicity 
and entropy dissipation were first reported in [22J for arbitrary N. However, computational 
studies have focused primarily on properties of the M\ model and its extensions, including 
multigroup equations and partial moment models (231124) . In related work, one may find 
simulations of M\ models based on other statistics, including Maxwell-Boltzmann (9Ullj and 
Fermi-Dirac [H[58]. This attachment to M\ is due to the fact that the higher order mem- 
bers of the Mfq hierarchy require the repeated solution of expensive numerical optimization 
problems. However, simulations of the M2 model [44(l63j (the next member in the hierarchy) 
have been performed for Bose-Einstein statistics and for Mjy up to order N = 15 for special 
benchmark problems using Maxwell-Boltzmann statistics (Tll28j. 

There are several reasons to consider perturbative modifications to EB models. First, it 
is more economical to improve the model with perturbative corrections than to increase the 
number of moments, since the latter increases the memory footprint and makes the defining 
optimization problem more difficult to solve. In this case of grey photon transport, the flux 
in the Mi model can be expressed analytically, i.e., no direct solution of the optimization 
problem is required. Thus, in this case, the argument against increasing the number of 
moments is especially compelling. A second reason is that perturbations add (among other 
things) diffusive terms to the EB model. It is hoped that these terms will smooth out non- 
physical shocks which are known to exist in EB models. These shocks are a generic artifact of 
the modeling procedure that result when approximating linear transport in phase space by a 
nonlinear hyperbolic balance law for a set of moments. A third reason is that the specification 
of boundary conditions for moment equations that are consistent with the underlying kinetic 
boundary conditions is an open problem. However, at least in the case of linear moment 
equations, recent efforts [37] have shown the potential for well-posed boundary conditions for 
models with perturbative corrections. A fourth and final reason is that entropy-based closure 
do not depend on the properties of the material. Perturbations on the other hand can couple 
material properties into the closure. 

Our goal in this work is to assess the qualitative behavior of the PM\ model relative to 
the original M\ model. We consider several test cases and find that the PM\ model gives 
mixed results. Roughly speaking, it does quite well for shock problems that the M\ model 
cannot handle. However, for more regular solutions, the two models perform comparably; 
and in some cases, the M\ model performs slightly better. 

One of the fundamental questions associated with any moment model is the issue of 
realizability. In the context of the M\ and PM\ models, the two unknowns E and F are called 
realizable if and only if they are the first two moments of an underlying kinetic distribution. 
This requirement on E and F is mathematically equivalent to the condition 



which must be satisfied point-wise in space and time. Here, c is the speed of light. It is 
expected that the solutions of the M\ model will satisfy (pQ) because it (like all EB models) 
is derived assuming an ansatz for the kinetic distribution which is positive. However, the 
underlying ansatz for the PM\ model is a perturbation of the EB ansatz that is no neces- 
sarily positive. Therefore, a modification of the PEB ansatz is needed which controls the 
contribution of the perturbative term. 

Even for the M\ model, the realizability condition ([1]) can be destroyed by a numerical 
method unless special care is taken to enforce it. To address this issue in the current setting, we 
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build on previous work with the M\ model [15] , using a Runge-Kutta discontinuous Galerkin 
(RKDG) method that is equipped with a special slope limiter [HS1EZ] in the spatial variable. 
For implementation of the PM\ model, this special limiter must be applied in combination 
with a control parameter that limits the size of the perturbations in the underlying ansatz 
of the PEB closure. The RKDG method [6] is a natural discretization here because we deal 
with a hyperbolic system of equations that is augmented by a diffusive term. 

The remainder of the paper is organized as follows. In Section[2j we introduce the radiative 
transfer equation and moment model framework. In Section [3l we derive perturbed entropy- 
based closures and give explicit expressions for the perturbed M\ model. In Section [5j we 
focus on the PM\ model in slab geometry and give details of the discontinuous Galerkin 
method used for simulation. In Section [61 we present numerical results. Section [7] is for 
discussion and conclusions. Several calculations and proofs are relegated to the Appendix. 



2 Radiative Transfer and Moment Equations 

We consider a collection of photons which move at the speed of light c through a static material 
medium. In engineering and physics applications, the fundamental quantity of interest is the 
radiation intensity tp = tp(x,Q,,v,t) which is a function of position x E K C M 3 , direction 
n, E § 2 , frequency v E (0, oo), and time t E (0, oo). Roughly speaking, tp is the flux of energy 
through a surface normal to Q. If / is the kinetic density of photons — that is, the number 
density with respect to the Lebesgue measure dxd^ldv — then tp = hvcf, where h is Planck's 
constant. 

The material is characterized by a temperature T = T(x), an equation of state for the 
energy e = e(T), and by scattering, absorption, and total cross-sections: a s , <7 a , and cr t = 
<r a + <r s that depend on x directly and also indirectly through the material temperature. 

2.1 The Radiative Transfer Equation 

The radiative transfer equation, which approximates the evolution of tp, is given by 

-d t tP + n-V x tP = C(t/j;T) . (2) 
c 

The collision operator C models interactions of photons with the medium. For the purposes 
of this paper, we assume C has the form 

C(V; T) := -a t i> + -3- {a s (j) + <r a £(T) + s) , (3) 

47T 



where <f> is the angular integral of tp: 



and the Planckian 



tpdn , (4) 



, , 2hv 3 1 

models blackbody radiation from the material. The constant k is Boltzmann's constant. The 
first term in C accounts for the loss of photons at a given frequency and angle due to both 
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out-scattering and absorption by the material. The second group of terms gives the gain of 
photons due to in-scattering from other angles, re-emission by the material, and a generic 
external source s. To make calculations, s is assumed to be isotropic. However, such an 
assumption is not necessary. 

The evolution of the material temperature is determined by a balance of emitting and 
absorbed photons: 

d t e(T) = a a ((V) - acT 4 ) , (6) 
where angle brackets are used as a shorthand notation for integration over angle and frequency: 



(.)= / / (-)dndu, (7) 

Jo Js 2 

and the T" 4 term in the first equation comes from the Stefan-Boltzmann Law: 

/ B(T)dv = acT 4 . (8) 
Jo 

The constant a is the radiation constant. Though the material equation ([6]) plays an important 
role, we will focus here on simulating the transport equation ([2|). 



2.2 Moment Equations 

The large phase space on which ([2]) is defined makes direct numerical simulation prohibitively 
expensive. Thus, approximate models are needed to reduce the size of the system. A common 
and well-known approach is the method of moments, for which the angular and/or frequency 
dependency of tp is approximated using a finite number of weighted averages. 

Derivation of any moment system begins with the choice of a vector-valued function 
m : S 2 — > M. n , Q, i— > [mo(O), . . . , m n -i(Q)] T , whose n components are linearly independent 
functions of £1. Evolution equations for the moments u(x,t) := (imfj(x, •, t)) are found by 
multiplying the transport equation by m and integrating over all angles to give 

-9 ( u + V I -M = (m^;T)). (9) 

c 

The system ([9]) is not closed; a recipe, or closure, must be prescribed to express unknown 
quantities in terms of the given moments. Often this is done via an approximation for ijj in 
([9|) that depends on u, 

ip{x,Q,t) ~ £(u(s,i))(«) , (10) 
and satisfies the consistency relation 

(m£(u))=u. (11) 

The resulting moment system is 

-dtu + V x ■ (flm£(u)} = (mC(£(u); T)) . (12) 

c 

In general, a closure is required to evaluate both the flux terms and the collision terms in 
(|9|). However for the collision operator in ([3]), no closure is required. Indeed, it is straight- 
forward to show that (mC(£ (u); T)) = (mC(^;T)) for any reconstruction that satisfies the 
consistency relation. Thus, for the purposes of this paper, we will be focused on closure of 
the flux term. As one might expect, the behavior of a moment system — and in particular 
its ability to capture fundamental features of the kinetic description — depends heavily on the 
form of the reconstruction. 
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3 Entropy-Based and Perturbed Entropy-Based (PEB) Clo- 
sures 



In this section, we briefly review the theory of entropy-based closures for radiative transfer 
} \63\ and introduce our new perturbative model. 



3.1 Entropy-Based Closures 

A general strategy for prescribing a closure is to use the solution of a constrained optimization 
problem 

min %(g) (13) 

s.t. (m§) = (rm/>) (14) 

where %{g) := (g(g)) and rj : R — > R is a strictly convex function that is related to the 
entropy of the system. For photons, the physically relevant entropy comes from Bose-Einstein 
statistics and is given by 



2ku 2 

Vid) = —3- K lo g( n ff) - ( n g + !) l °s(n g + 1)] , (15) 



c 

where the n g is the occupation number associated with g: 

„2 



n 9 := ^g . (16) 



The solution of (|13p is expressed in terms of the Legendre dual 



»/.(/) = — ^-logfl-expf-— /) ) . (17) 



Let 

B(a) := V i (a T m) = ^ -. . 1 - N . (18) 

V ' '* V ; c 2 exp(-^a T m) - 1 17 

Then we have the following. 

Theorem 1. The solution of (|13p is given by B{a), where a = a(u) solves the dual problem 

{(??* (a r m)> -a T u} . (19) 



mm 



It is also the Legendre dual variable of u with respect to the strictly convex entropy h(u) :- 
n{B{a{u))), i.e., 



ct(u) 



dh ^ 
ou 



(20) 



The moment system derived by setting £(u) = B(a) in (j 12 j) is hyperbolic and symmetric 
when expressed in the a variables and its solution formally dissipates h. Moreover, £ is an 
inherently positive quantity. 



Proof. The form of the minimizer in (|18|) can be derived formally using standard Lagrange 
multiplier techniques. However, a rigorous proof requires more technical arguments which 
can be found, for example, in [33] for the Maxwell-Boltzmann entropy and applied directly 
to the current setting. Once the existence of a minimizer is found, the other properties can 
be verified, as is done in [22l[35] □ 
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3.2 Perturbed Entropy-Based (PEB) Closures 

Perturbations to standard P/v closure^] have been derived for N = 3 in [37J and for general N 
in [53] (see also [27J and [60] ). The idea behind the derivation in [53] is to write tp = tp pn + ip, 
where tfj pn is the standard Pjy expansion. The perturbation ip satisfies its own kinetic equation, 
which can be then used to approximate tp in terms of ^ pn . The resulting "Dat" models gain 
a diffusive term in the equations for the highest order moments, uch an approach need not 
be restricted to the Pn equations. Indeed, following this exact strategy, we define 

1. The moment map Ai : g i— > u := (mg); 

2. The expansion map £ : u i— > £>(a(u)); 

3. The reconstruction 1Z = £ o M; 

4. The kinetic perturbation tp = tp — lZ(ip). 
The kinetic equation for ip is 

d t i> = d t i> - d t n(ip) = - d t s(u) = - s'(u)d t u 



where 



£'(u) = B'(a)^ = m T W(u) <mm r W(u) 



W(u) := rfl{a T m) 



and we have used the relation 



2h z v 



2,, 4 



exp(-^d T m) 



kc 



exp(- 



> 



Id = (m£») = (mB'(a)) |^ = <mm T W(u)> ^ 



9u 



to compute the matrix in ()22|) . By operating with "Pu := X — V u on 
£'(u).M, we can write (|21|) as 



1 



dt^ + V u (n-V x ^)=V u C(^T). 



It should be noted, for future use, that the projection Q u , given by 

1 



Qug 



■P u (W(u) 5 ) , 



(21) 

(22) 
(23) 

(24) 

where V u : = 
(25) 

(26) 



is self-adjoint in L 2 with respect to the positive weight W(u). 

Equation (|25p for the perturbation is exact. To derive a closure, we neglect the time 
derivative and perturbative component of the flux to arrive at the following approximate 
balance equation 

V u (fl-V x £(u)) ~P U C(V;T), (27) 

where 



V u C(ip;T) = -a t V u £(u)+^ 



1 

+ -r- 

47T 



<y s f^ + <y a P u B(T) + v u s 



(28) 



1 These closures are based on a spherical harmonic expansion in angle and can be formulated as an entropy- 
based closure with an L 2 cost functional I30U50I. 
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In the appendix, we show that for the grey equations with Bose Einstein entropy, 

P u £(u) dv = 0. (29) 



Knowing that this component will be integrated out in the final closure, we therefore solve 
(|27p for P u £(u) + tp in terms of a convective component ij) c and a diffusive component ■0 d : 



r u £{u) + v 



1 

47T 



r a P u (/> + r a P u B(T) + —P u s 



1 

0"t 



P u (n ■ V x £(u)) =: ^ + ^ d , (30) 



where r s and r a are the scattering and absorption ratios, respectively: 

r s = — and r a = — . 

Inserting (|30p back into the flux term of the moment equation (|12p gives 

(nmip) ~ (fimf (u)) + (OmV> c ) + (ftrm/> d ) =: f £ + f c + f D . 



(31) 



(32) 



At this point, it is not clear whether this flux dissipates an entropy or if the convective 
flux f c is always hyperbolic. In general, the hyperbolicity of moment models closed by an 
entropy minimization principle follows from the fact that (in terms the Lagrange multipliers 
a) the model can be written as a symmetric Lax-Friedrichs form [35] . This structure is not 
present here. However, at least for slab geometries, the convective flux in the PM\ model is 
hyperbolic. (See Proposition [2] in the following section.) Moreover, in general, the diffusive 
flux satisfies a local dissipation law. 

Proposition 1. The diffusion term f D dissipates the entropy h(u) := %(£(u)) locally in 
space. 



Proof. A dissipation law for h is found by multiplying the closed moment system (|12p by 
a T = Multiplying V x • f D on the right by a T gives 

a T (V x ■ f D ) = -a T [v* • (flma^Vu (ft • V x £(u))^ 

= -V x ■ (ni&^a^Vu (0 • V x £(u))\ + (v x a T ) ■ (Oma^V^ (ft ■ V x £(u)) 



where V x acts on the components of ft and the Lagrange multiplier a T on m. We only need to 
work with the term that is not in divergence form. We use the fact that B(a) = (hvc/k)m T yV 



(v a .a T ) • (sima^Pu (ft • V^(u))) = (v x a T "j • (sima^WQu 



hue 
~k~ 

hvc 



ft ■ v x £(u) 
w 



ft ■ V x (a T m)a- 1 WQ VL (ft • V x (a T m 



> 



where Q u := Id — Q u and Q u is given in (|26|) . 



□ 
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3.3 Controlling the Perturbations 



While the entropy-based ansatz in (|18p is positive for all $7, the addition of the perturbation 
in (|30p may lead to an ansatz which is not. As a consequence, the moments of the perturbed 
ansatz may not satisfy the realizability condition ([1]). To correct for this defect, we introduce 
a modification and approximate ip with 

S(u) = B(&) + 8^, (33) 

where 5(x,t) is a scalar control parameter. Several different choices for 5 are possible. For 
example, one could select it to ensure that S(u) is positive everywhere. However, this choice 
requires pointwise evaluations with respect to f2 — a task we would like to avoid. Instead, 
we select 5 in such a way as to preserve (pQ) in the numerical computation. While the exact 
form of 5 depends on the details of the numerical method, the general framework relies on 
the realizability conditions for the moments. We call an array (*S>o, Vl/i, . . . , ^n) realizable 
with respect to (1,0,..., Q® N ) if there exists a non-negative measure on d£ldv with density 
W(Q,v) such that *S> k = {Q® k *$>) for k = 1, . . . , N. The set TZn of all such vectors is called 
the realizable set. 

Roughly speaking, we select S to ensure that (^q, ^lj ■ ■ ■ ,^n) £ Note that such a 
5 always exists: When 5 = 0, there is no perturbative term and since the minimum entropy 
ansatz is always positive, (^o, . . . , \P/v) £ T^N- Details for the PM\ are given in Section 

E3J 



4 The Perturbed M x (PM X ) model 

The perturbed M\ model is based on the moments 

-(:)-(?)-($)■ 

where E is the photon energy density and F is the energy flux density. The model approxi- 
mates the evolution of E and F with the following system: 

d t E + V x • F = -aJycE - acT 4 ) + S, (35a) 
d t F + c 2 V x -IL(E,F) = -ca t F, (35b) 

where S(x, t) := J °° s(x, u, t) dv and the closure for the pressure term is 

U(E,F) := -((Ovn)(f(u)+f + ^ d )) =:U M ^{E,F) + U C (E,F)+U D (E,F) . (36) 
c 

Here n Ml (u) is the term that comes from the entropy ansatz (the entropy-based term). The 
term n c (u) is the convective correction and n D (u) is the diffusive correction. These correc- 
tions can be expressed in terms of n Ml and 

Q Ml := (O v3 £(u)> (37) 

which, in turn, can be expressed in terms of the unit vector n := F/\F\ and the scalars 

((O • n)*g> 

Xk ^ • (38) 
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Lemma 1. The correction terms II C and II D are given by 



n 



D 



ca t 



-V x -Q Ml + ^_(V x -F) + c ; 



n c = r) ■ [r s E + r a aT 4 + 



<9£ 
S_ 

ca t 



r)TT M i 

(v.-n^) 



<9F 



where n 



Proof, see appendix. 



(39a) 
(39b) 
□ 



Remark 1. T/ie formula for the convective correction is independent of the specific form of 
£ . In particular for the Pi model, the pressure term is n Pl = \E so that ^gi = lid and 
II C = 0/ cf. (|39bj) . This is consistent with the fact that the "D/v" models in \51$ contain only 
diffusive corrections. 

Lemma 2. The entropy-based terms II Ml and Q Ml are given by 



n Ml = |[(l-X2)Id+(3 X 2-l)(nVn)], 

Q Ml = ^[(Xi " Xs)(H V n) + (5 X3 - 3xi)n v3 ] 



where the scalars Xi> X2, an d X3 are defined in ([38 
Proof. See the appendix. 



(40a) 
(40b) 

□ 



In slab geometry, we end up with the following expressions for the components of the 
pressure term II = II Ml + II C + II D which can be computed from Lemma [1] and Lemma [2) 



n 



Mi 



X (E,F)E, n u = rj (E,F)[r s E + r a aT 4 + —) 

ca t J 



i 



— [D E (E,F)d x E + D F (E,F)d x F] =: D(u)cU 



ca t 



where the convection and diffusion coefficients are given by 

1 + 37 2 87 2 



(41) 
(42) 

(43) 
(44) 
(45) 

(46) 

2cE + ^/(2cE) 2 - 3F 2 

These coefficients are displayed in Figure [TJ Note that x> Vi an d Dp are all even functions 
of the ratio F/(cE), while Dp is odd. 



D E (E,F) 
D F (E,F) 



x(E, F) — 3 + 72 , 

3(7 2 + 5) ( 7 2 -l) 2 
2 7 4 ( 7 2 -3) 2 

9(7 2 + l)(7 2 -l) 2 
2 7 5( 7 2 _ 3 )2 



( 7 2 -3)ln 
( 7 2 -3)ln 



3(3 -7 2 )' 
1-7 



1 + 7 
1-7 
1 + 7 



67 
67 



and 



-3F 
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F/(cE) 



0.5 



-0.5 





F/(cE) 



0.5 



(a) Convection coefficients. 



(b) Diffusion coefficients. 



Figure 1: Perturbed Ml model coefficients. Left: x (dark green solid line) and 77 (black 
dash-dot line). Right: De (blue solid line) and Dp (red dashed line). 

5 Numerical Simulation using Discontinuous Galerkin 

In slab geometries, the diffusion-corrected M\ model and the material energy (|6j) reduce to 



d t u + d x f(u, d x u) 



d t T=^(E 



S(u), (x,t) e (X L ,X R ) X (0, tflnaj), 

aT 4 ), 



where 



u 



cE 
F 



s(u) 



-c 2 a a (E - aT 4 ) + cS 
-ca t F 



f(u,d x u) 



' cF 

c 2 il 



(47a) 
(47b) 



(47c) 



= n Ml + (5(II C + n D ) and C v = is the specific heat at constant volume. In this setting, 
the convective flux of the PM\ model is hyperbolic. 

Proposition 2. The perturbed M\ system in slab geometry is hyperbolic if I1 D = and 
\F\ < cE. 



Proof. The proof is a direct calculation, given in the appendix. 



□ 



We simulate the system (|47p using a Runge-Kutta discontinuous Galerkin (RKDG) method. 
The RKDG method is a method of lines: the DG discretization is only applied to spatial vari- 
ables while time discretization is achieved by explicit Runge-Kutta time integrators. The 
presentation here is rather brief and relies on details found in |48j . where the method was 
applied to the Mi model. A general description of the RKDG method can be found, for 
example, in [i~6|[T7]. 



5.1 Spatial Discretization 

We divide the computational domain (xl,xr) into J cells with edges 

XL = X l/2 < < . . . < Xj +1/2 = XR, 
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and let Xj denote the center of each cell Ij = (xj_i/2, Xj+x/2)- We let hj := Xj+i/a — x j-i/2 
be the length of the interval Ij and set h := m&Xjhj. We denote the finite-dimensional 
approximation space by 

V* ={»E L\x L ,x R ) : v u , E V k {I,), j = 1, . . . , J}, 

where V k (Ij) is the space of polynomials of degree at most k on the interval Ij. 

The semidiscrete DG scheme is derived from a weak formulation of (|47p . However, fol- 
lowing [18] we first reduce the convection-diffusion equations (|47p to a system of first-order 
equations by introducing the auxiliary variable v: 

d t u + d x f(u,v) = s(u), (48a) 
Sell = v, (48b) 

d t T= C -^{E-aT A ). (48c) 

The exact solutions u(-,i), v(-,i) and T(-,t) are then replaced by approximations u^(-,t), 
v/j(-, t) G x V h fc and Th(-,t) S V h fc , and the resulting set of equations is required to hold for 
all test functions bh £ V*: 



b h (x)d t u h (x,t)dx - f(u h (x,t), v h (x,t))d x b h (x)dx (49a) 
Ji 3 

+ {¥b h (x)j j = jf s(u & (z, t))b h (x)dx 
/ b h (x)v h (x,t)dx + / u h (x,t)d x b h (x)dx - lub h (x)j ■ = (49b) 
/ 6 /l (x)^T h ( 3 ;,t)dx= / ^(x)^^(^(x,t)-ar h 4 (x,t))^. (49c) 

J Ij Jlj V 

Here we use the bracket notation: 

[fb h (x)jj = £,-+1/26*^7+1/2) - { j-i/2h{xf_ l/2 ) (50) 

where 

f j±l/2(u,v) = f(u(^±i/ 2 ,t), v(x i±1/2 ,t)) (51) 

and 

6/1(^7+1/2) = lim , b h(x j+1/2 - e), b h (xf 1/2 ) = lim bhixj^li + e) (52) 

J+ ' e-5>0+ ' ' e^0+ 

are the right and left limits of bh at the cell interfaces 2^+1/2 ■ The term [116^(3;)] ■ is defined 
in an analogous fashion. 

Since the components of u/j(.,t) and v/ l (.,t) are piecewise polynomials, the edge values of 
u and v in (|5ip are not strictly defined. Thus, the nonlinear flux function f is replaced by a 
numerical flux f which depends on the pointwise limits of u^, on either side of the edge 
at Xj±i/ 2 : 

f j±i/2 = H u h(xj ±1/2 ,t),u h (x^ ±1/2 ,t),v h (xj ±1/2 ,t),v h (x+ ±1/2 ,t)). (53) 
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F/(cE) 



Figure 2: Eigenvalues of the hyperbolic flux Jacobian: Ml model (blue lines), perturbed Ml 
model (red lines). 



The notations for u carry over analogously. 

It remains to choose suitable numerical fluxes f and u. Since ()47|) has both a convective 
flux 

cF 

c 2 jjMj + F ) ( C(JaE + C(TaaT 4 + ^ 



f°(u) :-- 



and a diffusive flux 



f D (u,v) 





c 2 D(u) • v 



(54) 
(55) 



the choice is not obvious. Several approaches have been presented in literature 
In [6], the prescription for the diffusive term is given by 



f D 



u 



j±l/2 



fD ( U j±l/2' V j±l/2) + fD ( U j±l/2' V j±l/2 



U j±l/2 + U j±l/2 



(56) 

(57) 



Combining this term with the Lax-Friedrichs flux for f c (u) gives the following total numerical 
flux: 



f i±i/2 



l r 



f(u. ±1/2 , v . ±1/2 ) + f(u+ ±1/2> v+ ±1/2 ) - A(u 



/2 



1/2J 



J±l/2 



U 7±l/2) 



(58) 



where A is the largest magnitude of any eigenvalue of the Jacobian associated with f c . These 
eigenvalues, in general, depend on material properties, the temperature T and the source 
term S. In contrast to the Ml model, they are not bounded by the speed of light c. For 
example, neglecting the temperature and source the maximum value is approximately 9.12 c. 
We instead use the smaller value of A = c, which is the particle speed in the transport equation 
and is consistent with the application of the control parameter to enforce realizability (see 
Section l3.3p . Figure [2] shows the comparison of eigenvalues for the Ml and perturbed Ml 
modeling when c = 1, cr s = 1, <7 t = 3, T = = 5. 
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The DG solutions u^, and are expanded in terms of local basis functions {&/}f_ 
for V k (Ij) in each cell Ij\ 

1=0 1=0 1=0 

The standard choice of basis for V{Ij) is generated by Legendre polynomials Pi that are 
defined on the reference cell [—1,1]: 

bj(x) = P l ( 2{x ~ x ^ , je{l,...,J},le{0,...,k}, (59) 



and normalized so that 



1 2 



Pi(y)P m (y)dy = ——S l>m . (60) 
i 2m + 1 



With £j(y) := xj + yhj/2, this gives a formulation defined on the reference cell: 

~ J f(u{(Z 3 (y),t)y h ^(y),t))d y P m (y)dy (61a) 

+ fj+i/2 - (-i)"Vi/2 = ^ j\(ui^(y),t))P m (y)dy, 



Vyv^(t) + ^ nj(t) C,, m - u i+1/2 + (-l)^..^ = 0, (61b) 

2=0 

d t Ti(t) = ( f ^m(y)^K(£;(y))) dy (6ic) 



2m + l 



aCk + I Pm{y)Tl\^{ylt)a a {y)dy, 



2C V 
where 

Ci, m = J Pi(y)d y Pm(y)dy. (62) 

The remaining integrals are calculated by a quadrature rule. Note that (|61bp can be solved 
locally for vtn(t) in each cell ij, which can then be substituted back into (|61al) . 

5.2 Time Discretization: Explicit SSP Runge-Kutta Schemes 

The purpose of high-order, strong stability preserving (SSP) Runge-Kutta time integration 
methods is to achieve high-order accuracy in time while preserving desirable properties of the 
forward Euler method (for a review, see [25]). In this paper, we only use explicit schemes, 
which compute values of the unknowns at several intermediate stages. Each stage is a convex 
combination of forward Euler operators and this usually leads to modified CFL restrictions. 

The equations in (|6ip form a system of ODEs for the coefficients Um(t) and Tm{t). For 
all j G {1, . . . , J} and m £ {0, . . . , k}, we write this system in the abstract form: 

d t u J m (t) = £i >m (uh~~ \ . . . , u£~\ u£, . . . , u{, u J +1 , u{ +1 ), (63) 
d t Tl(t) = C? T (El ...,E{X,...,Ti). (64) 
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Here, £u,m and C? T m are the respective right-hand sides of the ODEs. 

Let {t n }n =0 be an equidistant partition of [0, tfi na i] and set At := t^ na \/N. Let A denote 
the application of a generic slope limiter, and let irym be the orthogonal projection onto 
the finite dimensional space W 71 . Note that A is applied at every Runge-Kutta stage. The 
algorithm for the optimal third-order SSP Runge-Kutta (SSPRK(3,3)) method reads as 
follows: 

• For all j E {1, . . . , J} and m E {0, . . . , k}, set u4° = A{7ry™(uo)}. 

• For all n E {0, . . . , N - 1}, j E {1, . . . , J}, and m E {0, . . . , A;}, 

1. Compute the intermediate stages 

u^)=A{u^ + At£i im K (0) )} 

u^ 2 ) = A j|u£* + + lAt£i jm (ui«)} (65) 

u^ 3) = A {^< n + + ^At^u^)} 

2. Set ujl = uj \ 

For the sake of completeness, we also state the optimal fourth order scheme SSPRK(5,4) [25]: 

= A{u[ n) +0.391752226571890 At ££ im (u& (n) )} 
u4 (2) = A{0.444370493651235 u^ n) + 0.555629506348765 u^ (1) 

+ 0.368410593050371 At £i m (u4 (1) )} 
u^ 3 ) = A{0.620101851488403 u4 (n) + 0.379898148511597 u^ (2) 

+ 0.251891774271694 At C J u>m (u4 (2) )} (66) 
u4 (4) = A{0. 178079954393132 u4 (n) + 0.821920045606868 u^ (3) 

+ 0.544974750228521 At C^ m (u^ 3) )} 

u { ^ +1) = A{0.517231671970585u4 (2) + 0.096059710526147 u^ (3) 

+ 0.386708617503269 u4 (4) + 0.063692468666290 At ££ jm (u4 (3) ) 
+ 0.226007483236906 At C{ m (u^ (4) )}. 

Note that SSPRK(3,3) permits a timestep of the same size as forward Euler, while the 
SSPRK(5,4) method is less restrictive, allowing for a time step that is 1.508 times larger 
the forward Euler scheme. 

5.3 Limiters 

As in |48j . two types of limiters are used. The first is standard; it is used to suppress spurious 
oscillations and maintain stability. There are many such limiters available. In this paper, we 
apply the moment limiter from |12| , which is a modification of the original limiter in [7] . This 
limiter is applied to the variables u, but not the auxiliary variables v or the temperature T. 
Additional details can be found in |48j . 



14 



5.3.1 Realizability-Preserving Limiter 

The second limiter is a realizability-preserving limiter which is needed to ensure that the cell 
averages of E and F satisfy the realizability condition ([1]) at each stage of the numerical 
computation. The limiter is based on the work from [67] and [66] and is very similar to what 
was done in [48] for the M\ model. The major difference here is the addition of the control 
parameter 5. 

An essential ingredient of the realizability limiter is the Gauss-Lobatto quadrature set 

{xj-i/2 =x},x 2 j ,..., xf' x ,xf = x j+ i /2 } C Ij, (67) 

where, for a spatial reconstruction of order k, M is the smallest integer such that 2M — 3 > 
2k + 1. This condition on M ensures accuracy of the scheme [15]. The weaker condition 
2M — 3 > k ensures that the quadrature integrates elements of the approximation space Vj? 
exactly. 

The realizability limiter is defined in order to ensure that u^Xj) E TZ 2 at each point x^ in 
the quadrature set. However, we enforce the convexity condition indirectly by requiring the 
positivity of the intermediate quantified! 

cE + F , „ cE-F 
Q := and R := . (68) 

The inverse transformation that maps (Q, R) t— > (cE, F) is given by 

E = an d F = Q - R. (69) 

c 

An additional limiter is also used to enforce the positivity of the temperature reconstruction 
at each quadrature point. 

We now proceed to define the limiters. Let u^' n = {cEjf 1 , F^ n ) and T^ n be the approx- 
imations of u and T in cell Ij at time t n , and let u^' n and T^' n denote the modifications of 
u^' n and T^' n that are generated by the limiting. We assume that the cell average of uj! , 
which we denote by u^ ,r \ is realizable, i.e., u^' n E TZ 2 . We also assume that the cell average 

of T^ n , which we denote by T^ n , is positive. Let Q J ^ n (x) and Rj^ n (x) be the approximations 
of Q and R, respectively, and define limited variables by 

Qf (*) = ej^Qir (x) + (1 - ^)Qf , (70a) 

(*) = + (1 - , (70b) 

f$ n {x) = OifT^x) + (1 - 4»m n , (70c) 



2 The meaning of all subsequent subscripts, superscripts and adornments of Q and R will be inherited from 
analogous definitions for E and F. 
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where 



e 3 / 



mm 



mm 



mm 



e/2 





O j,n 

t> mm 




-e/2 


K h ~ 


R j,n 
-"'mm 




n 


7fS> 
1 h - 


rpj,n ' 
min 



Q 



min 



R ],n 



min Ql' n (x^ 



min i?{' n (x'h 

-1 n.-r n J 



min Tl^ix 1 /) 
e=i,...,M 11 J 



(70d) 
(70e) 
(70f) 



The parameter e > is chosen to maintain numerical stability with finite precision arithmetic; 
its value should be small relative to the magnitude of the variables in a given problem. The 
components of u^' n are then defined using (|69p . They satisfy the following property which is 
a key ingredient for maintaining realizability in the RKDG scheme. 

Lemma 3 ( 08]). 7/u£ n G TZ 2 (respectively: T£" > 0^, then u£ n (x^) G K\ := TZ 2 + [e,0] T 



(respectively: T^' n (x^) > 0) f< 



or 



5.3.2 Setting the Control Parameter 

We now define the control parameter S, discussed in Section 
the following result. 



Our definition is guided by 



Lemma 4 ( [H]). In the one dimensional setting, a necessary condition for (cE, F, ells) G IZ3 
is that 

(ci) n 5 < e, 

(C2) |F± C n 5 | < cE±F. 

Rather than to require (cE, F, ells) G 7^-3, we choose 5 G [0,1] to ensure the weaker 
conditions (CI) and (C2). More specifically for any (cE,F) G 1Z 2 , we set 



5{E,F) 



U (E,F), U c (E,F)+U B (E,F)>0, 
U(£,F), U C (E,F)+U D (E,F) <0, 



(71a) 



where 

5q = min 



Mi 



E-U 

n c + n D 



mm 



2F + cE + cn Ml 2F + cE + clT vl1 



Mi 



c |n c + n D | 



c |n c + n D | 



(71b) 



Lemma 5. For all (cE, F) G K 2 , := U Ml + S[U C + n D ] satisfies (CI) - (C2). 



Proof. The assertion (cE,F) G 1Z 2 implies (cE, F, cIl Ml ) G IZ3. It follows then that for 
II D = 0, conditions (CI) - (C2) are trivially satisfied. It remains only to show the following 
inequalities: 

cU s < cE and cU 5 > 2F - cE and cU 5 > -2F - cE. (72) 

These relations are easily verified by applying the definition of 6 and using the fact that 
(cE, F, cII Ml ) G 7Z%. □ 
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With 5 given by (|7ip . one can show that cell averages of U/j remain realizable and that 
the cell average of remains positive in a forward Euler step. Let 

lf' n Wi';: *t* (73) 



Lemma 6. Assume that 2M — 3 > k and for each i = 1, . . . , M , that 

e K 2 , Tf n > (74) 
and IlJ" satisfies (CI) and (C2). Assume further that At satisfies the following conditions: 

(Al) At < min 



i,...,M \ca t p 



(A2) At < min 



wih 



i=l,...,M \ c(l + W£a t ^h) 

(A3) At < min J —. 1 . 

-e=i,.,M\a a/ ac(Tl'yf 

where h := mm,- hj. Then after a forward Euler time step, 

u£ n+1 G7e 2 and Tj n+1 > 0. (75) 

Proof. We refer the reader to the proof of Theorem 3 in [18] for the Mi model, which relies 
exactly on the conditions (A1)-(A3) and (C1)-(C2). The only difference is that (CI) and 
(C2) are assumed in Lemma El while in [35] they are naturally satisfied by the M\ model. □ 

Theorem 2. The Runge-Kutta discontinuous Galerkin scheme which combines 

1. the space discretization in (|6ip . 

2. the limiters in (J70]), 

3. the modified pressure lis in (03) with control parameter 5 given by (f7T]l , 
4- a strong-stability-preserving Runge Kutta time integrator, and 

5. a sufficiently accurate Gauss-Lobatto quadrature 



preserves the realizability of the moments in the sense of cell averages. In particular, if the 
time step conditions (Al)-(A2) in the statement of Lemma [7J hold and if u^' n € IZ2, then 

—j,n+l 



Proof. Application of the limiters in (|7U|) ensures that the conditions of Lemma [U] hold at each 
stage in the SSP-RK scheme. Each successive stage is an application of the forward Euler 
operator to the current stage with an appropriately modified time step. Thus, the conclusions 
of Lemma [6] apply at the next stage, including the final stage, which gives u^' n+1 . □ 
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6 Numerical Results 



In this section, we present several numerical computations in slab geometry for a choice 
of test cases that are common for the M\ model. The goal is to compare and contrast the 
perturbed M\ model with the M\ model and to point out benefits and drawbacks. Benchmark 
solutions are generated by the discrete ordinates method with an upwind scheme in space, 
high-order spherical harmonics or semi-analytic expressions. The RKDG implementation 
has been verified and benchmarked in [38]. The correct implementation of the additional 
perturbative terms has been checked by the method of manufactured solutions [53] . 

As in [IB] our algorithm is implemented in MATLAB, and Gauss-Lobatto quadrature on 
[-1,1] is used. Additionally, the Runge-Kutta time integration methods as well as parameters 
for the admissibility limiter are applied in the same way. In order to satisfy the conditions of 
Lemma [6] and to guarantee stability, the time step is set to 

At < min{ci,C2,C3,C4}, (76a) 

1 hw m i n C v h? 
c\ = , c 2 =-7— rr , c 3 = and c 4 = , (76b) 

CO- t ,max C(l + W ma xO"t,max/lJ aCT max 2(2fc + 1) 

where k is the polynomial degree, h = min^ hj, u> m in and u> max are the minimum and maximum 
quadrature weights, respectively. The quantities o"t jmax and r max are the maximum values of 
cr M and <7 a/ (T/' n ) 3 . 

The constant C4 is not needed to preserve realizability of the moments but rather to 
enforce a parabolic CFL condition. Without this condition, unstable modes grow without 
bound until the control parameter 5 turns on and damps them. Unfortunately, the parabolic 
CFL restriction leads to small time steps. 

The stability parameter for the realizability limiter of E and F is set to e = 10 -10 . The 
same value is also used to enforce conditions (CI) and (C2), i.e., the control parameter in 
rXTl is chosen such that 



cU s <cE-e and \F ± cILj| < cE ± F ± e. 

In Sections 16.1116.31 we study simulations with c = 1 and neglect the energy equation which is 
included in the last two cases from Section 16731 Unless otherwise stated, slope and realizability 
limiters are always turned on for all DG calculations. If transformation to characteristic 
variables for the slope limiter is used, it will be explicitly stated. 



6.1 Two-Beam Instability 

We consider two incoming beams at the boundaries of the domain [—0.5,0,5] and set c = 1, 
s = oE Particles stream from both boundaries in a purely absorbing material with a a = 4 = 
<7t and meet at x = 0. We avoid getting too close to the boundary of the realizability domain 
and represent these beams in our moment model using the boundary conditions 

u(0,t) = [1,0.9999] T , u(l,t) = [1,-0.9999] T , t>0 

and initial conditions 

u (x) = [2e,0] T , x G (-0.5,0,5). 

3 The radiation intensity of a beam is a delta distribution in angle at /i — 1 (left boundary) and \x = — 1 
(right boundary) which yields an energy density of E = 1/c and a flux density of F = ±1. 



IS 




-0.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 



x 



(a) t = 0.6 




-0.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 



x 

(b) t = 3 

Figure 3: Plots of E for the two-beam instability. J = 200, k = 2: M\ (purple circle 
perturbed M\ (blue dash-dot line), transport (black solid line). 
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For this problem, coupling to the material is ignored, and the material energy equation is not 
included in the simulation. 

In Figure El one can observe the formation of a shock in the Mi solution which persists 
at the steady-state. The perturbed M\ model also develops an unphysical transient profile 
in which the particle number jumps in the center of the domain. While this artifact persists, 
the steady state solution (t = 3) appears continuous. The steady-state solution also has 
noticeable kinks in the at x ~ ±0.3. For comparison, discrete ordinates solutions are plotted 
for which 256 discretization points in angle and 1000 points in space are used. The perturbed 
Mi is throughout closer to the transport solution. 

Remark 2. Precise explanations for the occurrence of shocks and kinks in the perturbed Mi 
solution require an additional analysis. For example, an explanation for the formation of 
shocks in the Mi model is given IlLtf . However, such an analysis goes beyond the purpose of 
this paper and will be postponed to future work. 

6.2 Source-Beam Problem 




0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 

x x 

(a) t = 0.5 (b) t = 1 




0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 



(c) t = 2 (d) t = 4 

Figure 4: Plots of E for the Source-beam problem. J = 300, k = 2: Mi (purple dashed line), 
perturbed Mi (blue dash-dot line), transport solution (black solid line) 
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In this problem, an incoming beam 

u(0,t) = [1,0.9999] T , t>0, 

on the left boundary of the domain [0, 3] hits an isotropic source S = 1/2 generating particles 
on the interval 1 < x < 1.5. In order to avoid complications of spatial discontinuities in 
the fluxes (HJEHI5S]; we smooth the source and material cross sections, which enter into the 
perturbative components of the flux. The source S is smoothed at the end points x = 1 and 
x = 1.5 



S 



i 

2' 

id 

0. 



PH 



1-A<x<l+A 
1 + A<x<1.5-A 
1.5- A<s< 1.5 + A 
else 



(77) 



Similarly, we design the material properties with the cross sections: 



<7-a 



1, 

(l-M^))/2, 
0, 

2, 

2 + 4(1+ P h(^)) 
10, 

10, 



< x < 2 - A 

2-A<x<2+A, 

else 

1-A<x<l+A 
1+A<x<2-A 
, 2-A<x<2+A 
2 + A < x < 3 
else. 



and 



(78) 



(79) 



The function pu is a Hermite polynomial of order 10 with p#(±l) = ±1 and (±1) = for 
k = 1,2,3,4. If ph is extended by ±1 respectively, it is a C 4 function. The material property 
functions are illustrated in Figure O 
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Figure 5: Source-beam problem. Material properties with A = 0.05: a s (red dashed line), a a 
(black solid line), S (blue solid line) 



On the right boundary, particles are absorbed and zero Dirichlet conditions 

u(3,t) = [e,0] T , £>0, 



(80) 
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are set. Initially, there are (almost) no particles in the system, i.e., 

u(x,0) = [e,0] T , s€(0,3). (81) 

The value of c is again set to one. 

The perturbed M\ results are compared to M\ and transport solutions. Classic M\ calcula- 
tions are performed using the DG method from [48] . with the same computational parameters 
as the PM\ model and slope limiting performed in the characteristic variables. The trans- 
port solution is computed using the discrete ordinates method with 600 spatial cells and 256 
discrete angles. 

One can observe in Figure 0] that as time increases, particles entering from the left bound- 
ary encounter the source in the interior. As this happens the M\ profile diverges from the 
transport solution. Even as steady state is achieved at t = 4, there is a large difference for 
x < 1. The PM± profile agrees much better with the transport solution. 



6.3 Gaussian Source 

The next test case simulates particles with an initial energy density that is a Gaussian dis- 
tribution in space and a zero energy density flux: 

T 

, £ = 0.1, xe(-L,L). 

Periodic boundary conditions on [— L,L] are prescribed where L = ifi na i + 1. The computa- 
tional domain is always chosen large enough to ensure that a negligible number of particles 
reaches the boundaries. No internal source is assumed (so 5 = 0), and the medium is purely 
scattering with er s = o"t = 1. The velocity c is also set to one and the material energy equation 
© is neglected. All DG results are computed with h = 0.01 and polynomial degree k = 2. 
For comparison, discrete ordinates solutions of the transport equations are obtained with 
h = 0.005 and 128 angular points. 

Figure [6] displays the solutions at ifi na i = 1,2,3,10. The Mi model gives the expected 
wave effects that are washed out at larger times. These effects do not occur in the perturbed 
Mi results. However, the PM\ forms Gaussian bell that are higher and more narrow than 
the benchmark solution. At lower times, their maximum propagation speed is roughly half 
the correct velocity. Nevertheless, at ig na i = 10 the front of the PM\ model catches up with 
the reference solution, at which point all three models agree reasonably well. 

The perturbation tp from Section ^. 2l is related to the difference between the Mi and trans- 
port solution. Figure [6] indicates that this quantity is highly time-dependent. Additionally, 
the spatial gradient of ip is large at shorter times. Hence, this numerical example violates 
the smallness assumptions made in the derivation of the perturbed Mi model in Section 13.21 
Thus the lack of accuracy is not surprising. 



u(s,0) 



£V2vr 



e W,Q 



6.4 Including the Material Energy Equation 

The next two examples involve ([2]) coupling to the energy equation ([6]). The linearized 
Marshak wave problem from [61] is analyzed first and then a Marshak wave with material 
parameters taken from |48j . 
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6.4.1 Smoothed Su-Olson's Benchmark Problem 

In [61], tabulated data is provided for analytic solutions to a linearized Marshak wave prob- 
lem, which serves as a validation of numerical algorithms in the radiative transfer community. 
In particular, this semi-analytic benchmark is also compared to diffusion-corrected Pjv ap- 
proximations in [53] • It is therefore of interest to study solutions of the perturbed M\ model 
to this problem. 

We compute approximations to ([2]) and (|6j) in slab geometry with the following physical 
data 

C„ = T 3 , c=l, <T a = l = 0- t . 

As in the source-beam problem, we seek to avoid the discontinuous material properties in [61] 
by introducing a smoothed version: 

' (1+Ph(^))/2, -0.5-A<x<-0.5 + A, 

I 1, -0.5 + A < x < 0.5- A, 

>- \ (!_ m (^5))/ 2) 0.5- A < x < 0.5 + A, 

0, else, 
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Figure 7: Plots of E for the Su-Olson Problem. J = 200, k = 2: Mi (purple dashed line), 
perturbed Mi (blue dash-dot line), semi-analytic (black asterisk). 



where pn is the Hermite polynomial described in Section 16. 21 However, we still compare our 
numerical results to the semi-analytic solutions from [61] because the length of the smoothing 
window 2A = 0.02 is relatively small. 

Initially, the medium is cold and there is no radiation: 

i/)(x,n,i/,t = 0) = and T(x,t = 0) = 0. 

Additionally, zero boundary conditions are enforced on an infinite domain: 

lim tp(x, fi, v,t) = and lim T(x,t) = 0. 

x— »±oo x— >±oo 

In practice, we impose periodic boundary conditions on a large domain [-L, L] where L = 

L^finalJ + 1- 

Solutions at different times are provided in Figure [7] for the half plane x > 0. A grid size 
of h = 0.01 and polynomial degree of k = 2 are chosen for all DG solutions. Classic Mi 
computations are slope limited in the characteristic variables. They are throughout close to 
the semi-analytic results. However at earlier times, the perturbed Mi solutions have larger 
slopes at x ~ 0.5. There, the perturbed Mi model yields larger deviations from the reference. 
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Only at t = 3.16228 solutions from both models are close to each other as well as to the 
semi-analytic points. 



6.4.2 Thin Marshak Wave 

In this problem, incoming radiation is prescribed on the left boundary by well-posed boundary 
conditions, 

T(0,t) = l, T(l,f) = 0, i>0, 

u(0, t) = [T(0, tfa, 0.8 • T(0, tfacf, u(l, t) = [2e, 0] T , t > 0, 
and the material is assumed to be purely absorbing: 

a a {T) = , <r s = 0, 5 = 0. 

(T + O.o)' 3 cm 

The physical constants are given by 

c = 3 • 10 10 cm/s speed of light, 

a = 1.372 • 10 14 erg/(cm 3 keV 4 ) radiation constant, 

C v = 3 • 10 15 erg/(cm 3 keV) heat capacity, 

which implies units of cm -1 for cross sections and keV for temperature T. Initially, the 
material is cold 

T (x) = 5- 10" 4 , s€(0,l), 
u(x,0) = [T (x) 4 a,0] T , x€(0,l). 

Due to above incoming radiation on the left boundary, radiation propagates through the 
medium from left to the right. The material temperature T (Figure [8^) and the energy density 
E (Figure [8^) decay smoothly to zero. The M\ and PM\ models yield very similar solutions. 
For comparison, a reference solution is computed using a P99 model that is calculated with 
the DG method from [12]. The simulation of the P99 model uses linear elements and 800 
spatial cells. The reference solution shows a much stronger decay in the energy and material 
temperature. 



7 Discussion and Conclusions 

In this paper, we have derived a hierarchy of closures based on perturbations of well-known 
entropy-based closures. The derivation has been done in the context of grey photon transport. 
Our derivations reveals final equations containing an additional convective and diffusive term 
which are added to the flux term of the standard closure. This is different to perturbations 
to standard P/v closures [54] which only gain a diffusive component. 

For the first member of the hierarchy, the PM\ model, we compute explicit formulas 
for all terms. The resulting system of equations is a convection-diffusion system which, for 
slab geometries, is discretized by using a Runge-Kutta discontinuous Galerkin method. By 
introducing a special limiter and an additional control parameter to modify the pressure term, 
we ensure that cell averages of the moments satisfy the important realizability property (pQ) . 
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(a) Temperature at t = 0.1ns. 

n 1 1 1 1 1 r 
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x 

(b) Scaled scalar flux a~ 1 E at t = 0.1ns. 

Figure 8: Thin Marshak wave. J = 160, k = 2: M\ (purple dashed line), perturbed M\ (blue 
dash-dot line), P99 (black solid circle line). 
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We perform simulations to compare qualitative results with the M\ model and with highly- 
resolved discretizations of the original transport equation. Improvements to the standard M\ 
model are observed in cases where unphysical shocks develop in the classical M\ model. 
However, for problems with continuous solutions, there is little or no improvement; and in 
some cases, the M\ model perform marginally better. 

Finally, we discuss some open problems in this framework which might be addressed in 
future: 

• Moment systems from entropy-based closures are known to be hyperbolic and satisfy 
a local dissipation law [22^35j. but such results are not yet known for the perturbed 
models. The partial result in Proposition Q] confirms that the additional diffusive term 
dissipates the entropy. Moreover, neglecting the diffusion contribution the system indeed 
becomes hyperbolic for the special case of the M\ model in one dimension (i.e., for slab 
geometry) . 

• Standard issues in analysis such as existence and uniqueness of solutions have yet to be 
investigated for the PM\ model or for PEB models in general. 

• In Section [5.3.21 the control parameter 5 is chosen to guarantee conditions (C1)-(C2). 
However, this ansatz is a crude modification of the original perturbative model and 
could distort numerical solutions. A more subtle approach, possibly along the lines of 
flux-limited diffusion, e.g. [39J, would be preferable. 

• An undesirable aspect of the RKDG method is the time step restriction required by the 
explicit time integrator. For convection-diffusion equations, stiff sources, and/or long 
time scales, this time step restriction is very harsh. To lower the computational effort, 
(semi-) implicit time discretizations are therefore necessary. In addition, the method does 
not address the challenges of spatially discontinuous fluxes that arise from discontinuities 
in sources and material cross-sections. 

• Another issue is the formation of unphysical shocks in the (perturbed) M\ solution. 
Further analysis is needed to understand why and when they appear and how to further 
mitigate them. 



A Appendix 

Computation of equation (|29p . We first calculate two frequency integrals. Let k := ^a T m 
and 9 := — kv. Then 



V* 







hue ^ \ . f°° 2hv A 1 



a 1 in | du = I — = ; : dv 



k J Jo c 2 exp(-Kv) — 1 



2" f" (82) 



c 2 k 4 J exp(<9) - 1 15c 2 k 4 
and 

' "(~ hvC a T m\ dv - [°° ex P(~ Kl/ ) dy 

V * \ k a m J V J kc [exp(-Kiz) - l] 2 

2 " 2 r *«pm m= ^X. (83) 



kcK 5 J [exp((9)-l] 2 15/cck 5 
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With these two integrals, it is easy to show that — a/4 is the unique solution to the linear 
system 



so that 



(mm T r]" ^— ^a T m^ \ @ = /m^ ^— ^OL T mj \ , (84) 

T „(-huc^ T / , ( -hvc „ T \\ a 

mm 77* — - — a m ) ( mr/ t — - — a m I ) = — — . (85) 



Using the definition of V u , 



4 

and again the integrals above: 



f°° -a T m 8vr 4 /i 2 -kn 8vr 4 /i 2 

J Pu (U) 4 15/fcc/c 5 " 4/ic 15fcc K 5 

Ah r°° 
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^(u) = ^r^a-m|, (86) 



15c 2 K 4 



y I — ~ j = y ^ 



□ 

Proof of Lemma The proof is a straight-forward calculation. It turns out to be more 
efficient to calculate I1 D and Il c without directly using (j22j) . Instead for any function g, we 
compute 

1 r)((Q\/0)F\ 1 f)T\ Ml 

<(£) V n)V u g) = ((ft V - £ U ) > (m k g) = <(0 V Q)<?> - ]T — (m fc5 ). (88) 

fc=o fc fe=o Ufc 

Using (j88|) . we find for the diffusive correction, 

n D = — ((nvnw 1 

ccr t \ 

= — 3- /(fi V 0)P U (0 • V x £ (u)) 



1 1 _* r)TT Ml 

— v x • <(o v3 £:(u)> + — V ^^v x . • ( mfc nf) 

fc=0 



fc=0 

1 r)TT M i 1 r)TT M i 

= -— v, ■ Q + — ^ (v. • f) + -^-(v, • n M - 

ca t ca t oE a t OF 

For the convection correction, we use the fact that <fi, B and S are independent of CI. This 
implies that (m^) = (Clef)) = and similarly for B and S. We also use the Stefan Boltzmann- 



28 



Law ([8]) and the identity uo = cE = J °° <\> dv. This gives 

c \ 

= it (<° v 0) M + £ ( (!! v + ( (!! v 0, ^ s 

r / r9TT Ml \ r f (9TT Ml 

= («" v - ^ir w ) + ss= («" v t! ' B » - ^r (B) 

i / <9n Mi 
5 u -^ L )(" 15 + "^ + ^)- < 90 » 

□ 

Proof of Lemma\M Let {ei,e2,e3} be any orthogonal basis for R 3 . Then 

3 3 

n = ^n iei , I2 i: =(n-ei), ^n? = i, (9i) 

i=i i=i 

and 

((Sl^u))^^^ £(11)^. (92) 

Now set e3 = n = F/\F\ and note that, according to Lemma [7] below, £(u) depends on 
only through f^. Thus, only the terms with even powers of £li and J7 2 will survive. For 
k = 2, this means 

C n Ml = (fii£(u)) eiVei + (^£(u)> e 2 V e 2 + (n|f(u)) n V n, (93) 

and for fe = 3, 

Q Ml = 3 (nlQ 3 £ (u)) ei V ei V n + 3 (n§n 3 £(u)) e 2 V e 2 V n + (J2|f (u)) n v3 . (94) 

The goal then is to write these formulas in terms of ^3 only. Let us focus first on IP* 1 . 
Because £(u) depends only on O3, symmetry arguments can be used to conclude that first 
two terms in (|93p are the same. Combined with the far right relation (|9ip . this gives 

C n Ml = (nl£(u)) (ei V ei + e 2 V e 2 ) + ((fi§£(u)) n V n 

= (nl£(u)) (ei V ei + e 2 V e 2 + n V n) + ((^3 - fif)£(u)) n V n 

= \ <(1 - n§)£(u)> Id + 1 ((30= - l)£(u)> n V n, (95) 

where we have used the fact that ei V ei + e 2 V e 2 + n V n is the identity. From the definition 
of \2i we conclude that 

n Ml = § [(1 - X2)H + (3 X2 - l)(n V n)] . (96) 
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(98) 



Similarly for k = 3, 

Q Ml = 3 (nln 3 £ (u)> (ei V ei + e 2 V e 2 + n V n) V n + ((«§ - 30?)0 3 £ (u)> n v3 

= | <(1 - 0|)0 3 £:(u)> Id V n + X - ((5^1 - 3)n 3 £(u)) n v3 
3rK 

= — [Ui - Xs)(Id V n) + (5 X3 - 3 X i)n v3 ] . (97) 

□ 

Lemma 7. For the M\ model, the multiplier 6t\ is co-linear with F, that is 

Wl\~W\ 

Proof. If £ (u) = rf* (— ^(ao + aimi)) solves the optimization problem (fl"3|) . then by defini- 
tion 

F=(nr)i (~(&o + \ • (99) 

Let R be any orthogonal 3x3 matrix which preserves F. Then multiplying (|99p by R gives 

Mlr£ (^-y(oo + «imi)jj) = (nri', (~^{&o + R&xiin)) \ , (100) 

where we have used the fact that the measure d£l is invariant under the action of R. Because 
the solution of the optimization is unique, we conclude that Ra± = a\ and therefore, since R 
is arbitrary, 6t\ and F must be co-linear hjhkjhkj □ 

Proof of Proposition [J| Without loss of generality, we consider c = 1 and prove that the 
eigenvalues of the Jacobian associated with the convective flux in (|54p are real. To do so, the 
following definitions are introduced: 

Q:= JL( {£ + ,, or < + |), e^e-fe+^ + l 

((f) := X(f) + r.vlf), f ~ F/E. 

We show that the radical a + /3 2 /4 in the formula for the eigenvalues is positive for all / ^ 1. 
Note that (|39bj) implies rj = 1/3 + x'f ~ X an d hence, 

t = r s (\-X + x'f)+X- (101) 



v 3 

The prime notation always refers to the derivative with respect to /. With this, we conclude 
/3 2 + 4a = e - 4/£' + 4£ = £' 2 - 4/£' + 4r s Q - X + x'f) + *X (102) 

= (£' - 2/) 2 + 4r s Q - X + X'/) + 4(x - / 2 )- (103) 
Using (|43p . straight-forward calculations imply 

X-/ 2 >0 for all //l and ^ - X + x'f > 0- (104) 
Applying (|104|) on (|103|) completes the proof. 

□ 
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